FRACTIONAL REACTION-DIFFUSION 
EQUATIONS 

O ■ R-K. SAXENA 

O . Department of Mathematics and Statistics, Jai Narain Vyas University 

^ Jodhpur, 342005, India 



(N 






m 



A.M. MATHAI 
Department of Mathematics and Statistics, McGill University 
Montreal, Canada H3A 2K6 
< 

U ; H.J. HAUBOLD 

^ ' Office for Outer Space Affairs, United Nations 

td '• P.O.Box 500, A1400 Vienna, Austria 



CN ' Abstract. In a series of papers, Saxena, Mathai, and Haubold (2002, 2004a, 



O 



2004b) derived solutions of a number of fractional kinetic equations in terms 



> 

t^^ , of generalized Mittag-Leffier functions which provide the extension of the 

^ ! work of Haubold and Mathai (1995, 2000). The subject of the present paper 

Q I is to investigate the solution of a fractional reaction-diffusion equation. The 

^ ■ results derived are of general nature and include the results reported earlier 



by many authors, notably by Jespersen, Metzler, and Fogedby (1999) for 
anomalous diffusion and del-Castillo-Negrete, Carreras, and Lynch (2003) for 
react ion- diffusion systems with Levy flights. The solution has been developed 
in terms of the H-function in a compact form with the help of Laplace and 
Fourier transforms. Most of the results obtained are in a form suitable for 
^ ■ numerical computation. 

1 Introduction 

Reaction-diffusion models have found numerous applications in pattern for- 
mation in biology, chemistry, and physics, see Murray (2003), Kuramoto 
(2003), Wilhelmsson and Lazzaro (2001), and Hundsdorfer and Verwer (2003). 
These systems indicate that diffusion can produce spontaneous formation of 
spatio-temporal patterns. For details, one can refer to the work of Nicolis and 



Prigogine (1977) and Haken (2004). A general model for reaction-diffusion 
systems is investigated by Henry and Wearne (2000, 2002) and Henry, Lang- 
lands, and Wearne (2005). 

The simplest reaction-diffusion models are of the form 

dN d'^N 

^ = ^^^ + ^(^).^ = ^(^.^). (1) 

where d is the diffusion coefficient and F{N) is a nonlinear function rep- 
resenting reaction kinetics. It is interesting to observe that for F{N) = 
'~fN{l — A^^), eq. (1) reduces to the Fisher-Kolmogorov equation and if we 
set F{N) = N[l — iV^), it gives rise to the real Ginsburg-Landau equa- 
tion. Del-Castillo-Negrete, Carreras, and Lynch (2002) studied the front 
propagation and segregation in a system of reaction-diffusion equations with 
cross-diffusion. Recently, del-Castillo-Negrete, Carreras, and Lynch (2003) 
discussed the dynamics in reaction-diffusion systems with non-Gaussian dif- 
fusion caused by asymmetric Levy flights and solved the following model 

dN 

— = vD:N + F{N), N = N{x, t), (2) 

with F = 0. 

In this paper we present a solution of a more general model of reaction- 
diffusion systems (2) in which ^ has been replaced by ^p^,/3 > 0. This 
new model extends the work of Jespersen, Metzler, and Fogedby (1999) and 
del-Castillo-Negrete, Carreras, and Lynch (2003). Most of the results are 
obtained in a compact form suitable for numerical computation. 

A generalization of the Mittag-Leffler function (Mittag-Leffler, 1903, 1905) 

oo n 

^-{z) ■■= E f? -TT' « e C, Re{a) > 0, (3) 

was introduced by Wiman (1905) in the generalized form 

oo n 

EaA^) ■= E w T^' a,(3eC, Reia) > 0. (4) 

The main results of these functions are available in the handbook of Erdelyi, 
Magnus, Oberhettinger, and Tricomi (1955, Section 18.1) and monographs 
by Dzherbashyan (1966, 1993). 



The H-function is defined by means of a Mellin-Barnes type integral in 
tlie following manner (Matliai and Saxena, 1978) 



where i 
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n™ r(6, + 5,0 iiU^{i-a,-A,i) 



i^j=m+l^\^ 



h,-B,i) m^^^,V{a,+A,i) 



(6) 



and an empty product is always interpreted as unity; m, n,p,q G A'^o with 
0<n<p,l<m<q, Ai, Bj G R+, Oj, bj G -R or C{i = 1, . . . ,p;j = 1, . . . ,q) 
such that 

Ai{bj + k) ^ Bj{ai - I - l),k,l e No;i = 1, . . . ,n;j = 1, . . . ,m, (7) 

where we employ the usual notations:A^o = (0, 1, 2 . . .), i? = (— oo, oo), R^ = 
(0, oo) and C being the complex number field. The contour Q is either 
L_oo, -^+00, or Li^oo- These contours are defined in the monographs by Prud- 
nikov, Brychkov, and Marichev (1989), Mathai (1993), and Kilbas and Saigo 
(2004). However, the explicit definitions of these contours are given below to 
complete the definition of the H-function. 

(i) Q = L_oo is a left loop situated in a horizontal strip starting at the 
point —00 + tip — 1 and terminating at the point —oo + iip2 with —00 < (^i < 

(P2 < +00; 

(ii) Q = -L+00 is a right loop situated in a horizontal strip starting at the 
point +00 + iipi and terminating at the point +00 + iip2 with —00 < ipi < 

(P2 < +00. 



I^iii) Q = Li^oo is a contour starting at the point 7 — ioo and terminating 



at the point 7 + ioo, where 7 G -R 



-00, +00 



A detailed and comprehensive account of the H-function is available from 
the monograph by Mathai and Saxena (1978), Prudnikov, Brychkov, and 



Marichev (1989), and Kilbas and Saigo (2004). The relation connecting 
p\l/g(z) and the H-function is given for the first time in the monograph by 
Mathai and Saxena (1978, p. 11, Eq.1.7.8) as 



\1> 



{ai,Ai),...,{ap,Ap) 
{bi,Bi),...,(bq,Bg) 



^p,q+l 



—Z 



(l-ai,yli),...,(l-ap,ylp) 
(0,l),{l-6i,Bi),.. .,(1-65,5,) 



where p\l/g(z) is Wright's generahzed hypergeometric function (Wright, 
1940); also see (Erdelyi, Magnus, Oberhettinger, and Tricomi, 1953, Section 
4.1), defined by means of the series representation in the form 



P%W 



V^<i 









(9) 



where z E C, ai, bj G C, Ai, Bj E R = (—00, 00), Ai, Bj j^ {i = 1, . . . ,p;j = 
1, . . . , g), J2'j=i Bj — J2^=i Aj > —1;C being the set of complex numbers and 
T{z) is Euler's gamma function. This function includes many special func- 
tions. It is interesting to observe that for Ai = B 
reduces to a generalized hypergeometric function pFq{z) as 

_ n,'=ir(a, 






(ap,l) 



n?=ir(6, 



J — l,Vi and j, eq. (9) 

<? 

pFg{ai,...,ap;bi,...,bg;z), (10) 



where aj 7^ —i^{j = 1, • • • ,P and u = 0,1,2, ... .);p < q or p = q,\z\ < 1. 
Prior to (9), Wright (1933) introduced a special case of (9) in the form 



^(a,b:z) 



0^1 



(6,a) 



E 

r=0 



T{ar + b) (r)! 



which widely occurs in problems of fractional diffusion. It has been shown 
by Saxena, Mathai, and Haubold (2004b) that 



EaA^) 
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l),(l-/3,a). 



If we further take /? = 1 in (11) and (12), we find that 
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(11) 

(12) 

(13) 
(14) 



where Re{a) > 0,a E C. 

From Mathai and Saxena (1978) and Prudnikov, Brychkov, and Marichev 
(1989, p. 355, Eq.2.25.3), it follows that the Laplace transform of the H- 
function is given by 



L {t^-'//™," 



zr 
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(l-p,cr),(ap,Ap) 



(15) 



where a > 0,Re{s) > 0, i?e[p + a]"<7<„(;^)] > 0,\argz\ < [n/2]9,e > 0; 

By virtue of the cancellation law for the H-function (Mathai and Saxena, 
1978), it can be readily seen that 
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(16) 



where a > 0, Re{s) > 0, Re[p + a]^ 



/l-O, 

<j<n\ Ai 



)] > 0, \argz\ < \'k6i,6i > 0; 



9 = 9 — a. Two interesting special cases of (16) are worth mentioning. If we 
employ the identity (Mathai and Saxena, 1978) 



1,0 



-"0,1 



X 



i»A) 



x°'exp{—x) 



we obtain 



L-^[s-Cexp{-zs'')] = tP-^Hli 



zt-"" 



(0,1) 



(17) 
(18) 



where Re{s) > 0, (T > 0. 

Further if we use the identity (Mathai and Saxena, 1978) 



R< 



2,0 
0,2 



X 



(f.l){-|.l) 



eq.(16) yields 



2K^{2x 



2+-20- 



1/2^ 



2L-\s-PK^{zs'')] =t''-^Hl^ 
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{P,2<7) 

{f,i),(-f,i) 



(19) 



(20) 



where Re{p) > 0,Re{z'^) > 0,Re{s) > 0, and Ki,{.) is the Bessel function of 
the third kind . 

In view of the result of Saxena, Mathai, and Haubold (2004a, p. 49), 
also see Prudnikov, Brychkov, and Marichev (1989, p. 355, eq.(2.25.3.2)), the 



cosine transform of the H-function is given by 



tP-^cos{kt)H^f' 



at" 






dt 



(21) 



T7-n+l,m 



A;^ 



(l-b„B,),{i±£,f) 
{p,^),(l-ap,A^),(i±£,$) 



where Re\p + Aii<j<m] > !> k^fi' '^1 < 2 ■^^^^S ^ > 0, 6* is defined with the result 
eq. (15). 

The Riemann-Liouville fractional integral of order v is defined by (Miller 
and Ross, 1993, p.45) 



1 (y) Jo 



(22) 



where Re{u) > 0. 

Following Saniko, Kilbas, and Marichev (1990, p. 37), we define the frac- 
tional derivative for a > in the form 



^Dtf{t) 



1 rf" /•* f{u)du 



V(n- a) df' Jo (t - m)"-"+i ' 



n= [a] + 1, 



(23) 



where [a] means the integral part of the number a. 
In particular, if < a < 1, 



oA"/W 



f{u)du 



dtV{l -a) Jo {t-uY' 
and if q; = n G A^ = {1, 2, . . . , } , then 

,D'lf{t) = D-f{t){D = d/dt), 



(24) 



(25) 



is the usual derivative of order n. 

From Erdelyi, Magnus, Oberhettinger, and Tricomi (1954b, p. 182), we 
have 

L{,Di^f{t)} = s-^F{s), (26) 

where 

F{s) = L {fit); s} = ris) = / exp{-st)f{t)dt, Re{s) > 0. (27) 

Jo 



The Laplace transform of the fractional derivative is given by Oldham and 
Spanier (1974, p.l34, eq.(8.1.3)) 

n 

L{,Dtf{t)} = s-F{s) -Y.S''' oDrfm=o- (28) 

r=l 

In certain boundary-value problems, the following fractional derivative of 
order a > is introduced by Caputo (1969) in the form 

m — 1 < a < m, Re{a) > 0,m E N. 

d'^f 
= — — , if a = m. (30) 

Caputo (1969) has given the Laplace transform of this derivative as 

m—l 

L {D^f{t); s} = s'^Fis) - ^ s"-"-V"(0+), m - 1< a < m. (31) 

r=0 

The above formula is very useful in deriving the solution of differ-integral 
equations of fractional order governing certain physical problems of reaction 
and diffusion. We also need the Weyl fractional operator defined by 



1 rf" /•* f{u)du 
r(n — fi) dt^ j-oo (t — u) 



-ooD^fit) - _ _ ^^_„+l , (32) 



where n = [fj] is an integral part of ;U > 0. 

Its Fourier transform is (Metzler and Klafter, 2000, p. 59, A. 11, 2004) 

F{^^D^J{x)} = {tkrf{k), (33) 

where we define the Fourier transform as 

h{q) = / h{x)exp{iqx)dx. (34) 

Following the convention initiated by Compte (1996), we suppress the imag- 
inary unit in Fourier space by adopting the slightly modified form of above 
result in our investigations (Metzler and Klafter, 2000, p. 59, A. 12, 2004) 

F{_^D^J{x)} = -\krm (35) 



instead of (33). Finally we also need the following property of the H-function 
(Mathai and Saxena, 1978) 
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x' 






_ TTm,n 
^ P,<1 
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(36) 



where (5 > 0. 



2 The Fractional React ion- Diffusion Equation 

In this section we will investigate the solution of the reaction-diffusion equa- 
tion (37). The result is given in the form of the following theorem. 
Theorem. Consider the following fractional reaction-diffusion model 

^|j'^^ = V ^ooD^N{x, t) + ^{x, ty,r],t>0,xeR,0<f3<2, (37) 

with the initial condition 

N{x, 0) = f{x), Nt{x, 0) = g{x) for x E R, (38) 

where Nt{x,0) means the first partial derivative of N{x,t) with respect to 
ip evaluated at t = 0,?7 is a diffusion constant and ip{x,t) is a nonlinear 
function belonging to the area of reaction-diffusion. Then for the solution of 
(37), subject to the initial conditions (38), there holds the formula 



N{x,t) = — / f{k)Ep^i{-r]\k\"t^)exp{-ikx)dk (39) 

zvr j-oo 



1 r°° 
+ — / tg{k)Ep^2iv\k\''t'^)exp{-ikx)dk 

+ Trh'''' ^{k,t - OEp^^{-v\kn^)exp{-zkx)dkdC 
zvr JO J-oo 



Proof. If we apply the Laplace transform with respect to the time variable 
t and Fourier transform with respect to space variable x and use the initial 
conditions and (32), then the given equation transforms into the form 

N*(k A_i(^>!li g{k)s^-^ ^*{k) 



On taking the inverse Laplace transform of (40) and applying the result 



^a-l 



-1 



t'3-"E^,^_,+i(-at^), (41) 



where Re{s) > 0, Re{P — a + 1) > 0, it is seen that 

N{k,t) = f{k)Ep,,{-r]\krtf') + ~g{k)tEp,2{-v\krt(') (42) 

+ f0{k,t-Of''EfsA-v\kn^)dt 
Jo 

The required solution (39) is now obtained by taking the inverse Fourier 
transform of (38). Thus, we have 

N{x,t) = — / f{k)Ep^i{-i]\k\''t'^)exp{-ikx)dk 
zvr J-oo 

1 /■°° 
+ — / tg{k)Ef3,2i-v\k\''t'^)ewi-ikx)dk 

271 Joo 

1 rt roo 

+ Trh"'' v*{k,t-OE(,A-v\kn>M-^kx)dkd^. 

In JO J-oo 



This completes the proof of the theorem. If we set a = 2, we obtain the 
result given by Debnath (2003). 

Note. By virtue of the identity (12), the solution (39) can be expressed in 
terms of the H-function as can be seen from the solutions given in the special 
cases of the theorem in the next section. 

3 Special Cases 

When g{x) = 0, then applying the convolution theorem of the Fourier trans- 
form to the solution (35), the theorem yields 
Corollary 1.1. The solution of fractional reaction-diffusion equation 

Q^N{x, t) - r]_^D2N{x, t) = ^x, t), x e R,t > 0,r] > 0, (43) 

subject to the initial conditions 

N{x, 0) = /(x), Nt{x, 0) = for a; G i?, 1 < /3 < 2, (44) 



where rj is a diffusion constant and f{x,t) is a nonlinear function belonging 
to the area of react ion- diffusion, is given by 

/oo 
Gi{x-r,t)f{r)dT (45) 

-oo 

+ fit -if-' rG2{x-r,t-0^{T,0drd^, 
Jo Jo 

1 f^^ 

Gi{x,t) = — / exp{-ikx)Ep^i{-7]\k\'^t'^)dk (46) 

zvr j-oo 



where 
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dk 
, q; > 0, 
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,q; > 0. 



(47) 



If we set f{x) = 6{x), v? = 0, g{x) = 0, where 6{x) is the Dirac-delta function, 

then we arrive at the following 

Corollary 1.2. Consider the following reaction-diffusion model 



d'^Njx, t) 



V 



,D'^N{x, t),r]>0,xeR,0<P<l, 



(48) 



with the initial condition N{x,t = 0) = S{x), where r^ is a diffusion constant 
and 6{x) is the Dirac-delta function. Then the solution of (44) is given by 



N{x,t) 
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(49) 



In the case (3 = 1, then in view of the cancellation law for the H-function 
(Mathai and Saxena, 1978), (49) gives rise to the following result given by 



Jespersen, Metzler, and Fogedby (1999) and recently by Del-Castillo-Negrete, 
Carreras, and Lynch (2003) in an entirely different form. 



For the solution of fractional reaction-diffusion equation 
d 



dt 



with initial condition 
there holds the relation 



N{x,t)=r]^^D^N{x,t), 



N{x,t = 0) =S{x), 
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(50) 
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(52) 



where a > 0. 

It may be noted that (52) is a closed-form representation of a Levy stable 
law, see Metzler and Klafter (2000, 2004). It is interesting to note that as 
q; — > 2, the classical Gaussian solution is recovered as 
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It is useful to study the solution (49) due to its occurrence in certain 
fractional and diffusion models. Now we proceed to find the fractional order 
moments of (49). 

Here we remark that applying Fourier transform with respect to x in (48) it 
is found that 

^iV(A;,t) = -r^|A;riV(A;,t), 

which is the generalized Fourier transformed diffusion equation, since for 
a = 2 and for /5 = 1, it reduces to Fourier transformed diffusion equation 



dN{k,t) 
di 



v\k\'N{k,t), 



being a diffusion equation, for a fixed wave number k (Metzler and Klafter, 
2000, 2004). Here N{x,t) is the Fourier transform witli respect to x of 

N{x,t). 



4 Fractional Order Moments 

We now calculate the fractional order moments defined by 

roo 

< |x(t)|'^>= / \x\^N{x,t)dx. 

Jo 

Using the definition of the Mellin transform 

M{q(t);s}= / t'-\{t)dt, 
Jo 

we find from (49) that 

< |a;(t)|'^ >= / \x\^N{x,t)dx 
2 /•«= 



< kWr >= 



x'-'H''' 



a Jo 
Applying the integral formula 
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dx = a"«0(-O, 



(55) 
(56) 

(57) 
(58) 

(59) 



where i<j<mR^i.^) < R^iO < '"^aa;i<j<„-Re(^^), \arg a\ < ^7i9,9 > 0, 6^ is 
defined in (6) and 0(— i^i the definition of the H-function (6). We see that 

(60) 



, . ^,^ 2, .,(5 r(-^)r(i + 5)r(i + ^ 



where Re{6) > —1 and Re{6 + a) > 0. 

Two interesting special cases of (60) are given below. 

(i) 

1 



As 5^0, then using the result 



rf^) 



z for z « 1, 
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we find that 



h<kwi'>=i. 



(ii) When a = 2, the hnear time dependence is 



5^2,a^2 < F('^)I 



2?7t^ 



r(i + /3)' 
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(63) 



5 Behavior of the Solution in Equation (49) 



Eq. (49) can be written in terms of a Melhn-Barnes type integral as 
1 1 



N{x,t) 



s/3^ 



a\x\ 2-Ku Jl r(i - f )r(i - f )r(s/2 
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jjlafP/c 



ds. 



(64) 



Let us assume that the poles of the integrand are simple. Now evaluating 
the sum of residues in ascending powers of \x\ by calculating the residues at 
the poles of r(l — s) at the points s = 1 + u {u = 0,1,2, . . .) and r(l — -) at 
the points s = a{l + u) (z/ = 0, 1, 2, . . .), it is found that the series expansion 
of the general solution (45) is given by 
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where < Re{v) < a, | ^ J- < 1. From (65), we infer that 



N{x,t) -^ A + B\x 



o-l 



clS X 



0, 



(66) 



where A and B are numerical constants. Further from the series expansion 
(61), it can be seen that the initial behavior is given by 



N{x,t) ~ ^\" ,, , "^ ., , forl<a<2. 



rfi^)bi"-i 



^i/2^2"t'3r(i - /?)r(f ; 



-, for < a < 1, 



(67) 
(68) 



with {(^^1 « 1. 

Further, if we calculate the residues at the poles of T{s/a) at the points 
s = —au {u = 0, 1, 2, . . .), it gives 



for I a; I > 1. From (69), it can be readily seen that 



r]tl^ 



(69) 



Nix,t) ~ —r for large |a;|. (70) 

\x\ 

6 Conclusions 

React ion- diffusion equations are modeling tools for the dynamics presented 
by a competition between two or more species, activators and inhibitors or 
production and destruction, that diffuse in a physical medium. When the dif- 
fusion, a homogenizing process, presents characteristic length or time scales 
that are different for the species, a morphological instability may trigger 
spatio-temporal pattern formation. This mechanism and its various versions 
is the Turing instability (Nicolis and Prigogine, 1977; Haken, 2004). The 
complexity of such equations and systems of equations limits the possibilities 
to derive analytic closed-form representations for their solutions (Wilhelms- 
son and Lazzaro, 2001; Kulsrud, 2005; Hundsdorfer and Verwer 2003). 

This paper attempts to derive analytic closed-form representations of 
solutions for a general class of fractional react ion- diffusion equations (eqs. 
(37), (43), (48)), which also provide as special cases respective solutions 
of standard react ion- diffusion equations. For this purpose the paper sum- 
marizes specific techniques for Laplace, Fourier, and Mellin transforms, re- 
sults for Mittag-Lefiier and Fox functions, as well as the applications of 
Riemann-Liouville, Weyl, and Caputo fractional calculus for tackling frac- 
tional reaction-diffusion equations. Closed-form solutions of the equations 
are given in terms of Fox's function and their behavior for small and large 
values of the respective parameter is derived. Fractional order moments and 
behavior of the fundamental solution as given in in eq. (49) are provided in 
detail. 



The methods used in this paper may also be apphed to nonhnear Fokker- 
Planck equations, taking into account the different physical meaning of reac- 
tion, drift, and diffusion coefficients, respectively (Tsallis, 2004; Tsallis and 
Bukman, 1996). 
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